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Dynamics of ideal fluid with free surface can be effectively solved by perturbing the Hamiltonian 
in weak nonlinearity limit. However it is shown that perturbation theory, which includes third 
and fourth order terms in the Hamiltonian, results in the ill-posed equations because of short 
wavelength instability. To fix that problem we introduce the canonical Hamiltonian transform from 
original physical variables to new variables for which instability is absent. 
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' The Euler equations describing dynamics of ideal fluid with free surface is a Hamiltonian system, which is especially 
simple if the fluid motion is potential, v = V<&, where v is the fluid's velocity and \& is the velocity potential. In this 
case jH |H [3 the Euler equations can be presented in the form: 

5?! dji _6H_ _SH_ 

PLh I dt ~ <5* ' 8t ~ Sn ' ^ ' 

Here z = r)(r) is the shape of surface, z is vertical coordinate and r = (x,y) are horizontal coordinates, $ = $| z _ 
is the velocity potential on the surface. The Hamiltonian H coincides with the total (potential and kinetic) energy 
of fluid. The Hamiltonian cannot be expressed in a closed form as a function of surface variables rj, but it can be 
presented by the infinite series in powers of surface steepness |V?/|: 

: H = H + H 1 +H 2 + .... (2) 

in 

Here Hq, Hi, H 2 are quadratic, cubic and quartic terms, respectivel y. E quations Q), are widely used now for 
numerical simulation of the fluid dynamics 0, la, @> S & S OH EH LUJ LLJ- Hj ■ These simulations are performing by 



the use of the spectral code, at the moment a typical grid is 512 x 512 harmonics. 

Canonical variables are used also for analytical study of the surface dynamics in the limit of small steepness. It 
was shown 0, 0, 0] that the simplest truncation of the series @ , namely 

m ; 

H = H + H U (3) 

, leads to completely integrable model - complex Hopf equation. In framework of this approach one can develop the 
self-consistent theory of singularity formation in absence of gravity and capillarity for two dimensions (one vertical 

. coordinate z and one horizontal coordinate x). 

However, use of canonical variables 77, has a weak point, which becomes clear, if we concentrate our attention on 
the complex Hopf equation, 

dt 2\dx) ' { ' 

which comes from Eqs. Q, EJ- Here 

* = Re{n) (5) 

and n is the analytic function of the complex variable x in a strip —h < Im(x) < 0, h is the depth of the fluid. The 
weak point is that Eq. (jJJ is ill-posed. A general complex solution of this equation is unstable with respect to grow of 
small short-wave perturbations. The same statement is correct with respect to more exact fourth order Hamiltonian 



H — Hq + Hi + H 2 , 



(0) 
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which is used in most numerical experiments. These experiments are easy becomes unstable: to arrest instability one 
should include into equations strong artificial damping at high wave numbers. Even in presence of such damping one 
can simulate only waves of a relatively small steepness (not more than 0.15). 

In this Article we show that these difficulties can be fixed by a proper canonical transformation to another canonical 
variables. It is remarkable, but the property of nonlinear wave equation to be well- or ill-posed is not invariant with 
respect to choice of the variables. 

In the present Article we demonstrate that there are new canonical variables such that the Eqs. (JTJ, (JBJ are well- 
posed if we consider the nonlinearity up to the fourth order in the Hamiltonian. We call these variables "optimal 
canonical variables" . We demonstrate in the present Article that the choice of the optimal canonical variables is unique 
provided we additionally require the Hamiltonian system to be free of short wavelength instability for largest possible 
steepness of the surface, i.e. for the largest possible nonlinearity. We conjecture that optimal canonical variables allow 
simulation with higher steepness compare with standard variables <J/, 77. We can also formulate a conjecture that the 
optimal canonical variables exist in all orders of nonlinearity. 



BASIC EQUATIONS AND THE HAMILTONIAN FORMALIZM 

Consider the dynamics of an incompressible ideal fluid with free surface and constant depth. Fluid occupies the 
region 

-h<z<r}(r), r=(x,y), (7) 

where (x, y) are the horizontal coordinates and z is the vertical coordinate. 
Viscosity is assumed to be absent and the fluid's velocity v is potential one, 

v = V$, (8) 

where $ is the velocity potential. Incompressibility condition, 

V-v = 0, (9) 

results in the Laplace Eq. 

A$ = 0. (10) 
The potential $ satisfies also the Bernoulli equation: 

*t + ^(V$) 2 +P + 5^ = 0, (11) 

where p is the pressure, g is the acceleration of gravity, and we set density of fluid to unity. 

There are two types of boundary conditions at free surface for Eqs. 1)10(1. (|ll(l . First is the kinematic boundary 
condition 



dr] 
~dt 



$ z - V77 • V$ 



v n Vi + (yW, (12) 



where v n = n • V<i> is the normal component of fluid's velocity at free surface, and n = (— V77, 1) [l + (Vr/) 2 ] is 
the interface normal vector. 

Second is the dynamic boundary condition at free surface 

p\ =aV-^=5==, (13) 



where a is the surface tension coefficient which determines the jump of the pressure at free surface from zero value 
outside of the fluid to p\ value inside fluid. 

1 I Z — Tj 

Boundary condition at the bottom is 



(14) 
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Eqs. (|10|l — <|14|) form a closed set of equations to determine the dynamics of free surface. 
The total energy, H, of the fluid consists of the kinetic energy T, and the potential energy, U: 

H = T + U, 



1 i / dr I (V<I>n/:. 



U=\g fr, 2 dr + a f [^l+(Vr?) 2 -l 



dr. 



It is convenient to introduce the value of the velocity potential at interface as 



(15) 
(16) 

(17) 



(18) 



It was shown in Ref. j^] that the free surface problem (|10|l — 114|l can be written in the Hamiltonian form Q), with 
the Hamiltonian H defined in i|15|) . 
Fourier transform, 



* k = -!- / exp(-ik • r)*(r)dr, 

2lT J 



is the canonical transform which conserves the Hamiltonian structure and Eqs. Q take the following form: 



dt S-f-u 



dt 



(19) 



(20) 



WEAK NONLINEARITY 



If a typical slope of free surface is small, |V?7| <C 1, the Hamiltonian H can be series expanded (see Eq. ©) in 
powers of steepness |Vry| which gives 0,0] : 



\ J [A k \V k \ 2 + B k \ Vk \ 2 }dk, 
A k = fctanh(/j/i), B k =g + ak 2 , jfc=|k|, 
Hl = ^J i k 1 ) ,k 2 *k 1 *k 2 '7k3'5(ki + k 2 + k 3 )dk 1 dk 2 dk 3 



Ho 



1 



2(2tt) 2 

where matrix elements are given by 



r(2) 
ki ,k2 ,k3,k 



^Sk^ks,^*^^ - J(ki • k 2 )(k 3 • k 4 )r/ kl r/ k2 
x?7ic 3 ?7k4'5(ki + k 2 + k 3 + k i )dkidk 2 dk 3 dk4 : 



4 = jAxAa (a 1+3 + A 2+3 + A 1+4 + A 2+4 



-^{k\A 2 + klAx), A 3 = A k] , Aj+i = A kj+kl . 
The corresponding dynamical equations follow from JU, 10, ^ljl . (I22|) . Ij23(l : 



— = -grj + ctAj/ + - [(A*) 



-(A*)(i*)j?- 



(i*)i[r?(i*)] 
Vr/(Vr] ■ Vr?) 



<3?j 
^7 



i* - V • [(W)//] - + i/ryi [77!*] I 



~A[7, 2 i*] +ii[ry 2 A*] 



(21) 
(22) 

(23) 



(24) 



(25) 
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where A is the linear integral operator which corresponds to multiplication on ktanh(kh) in Fourier space. For two 
dimensional flow, ty(x,y) — ^(x), rj{x,y) — T]{x), this operator is given by 

A = -^R (26) 

i r + °° f(x') 

Rf{x) = Tk RV - L sinh[(^ -i/( 2 .)] dX ' (27) 

where P.V. means Cauchy principal value of integral. In the limiting case of infinitely deep water, h — > oo, the operator 
A turns into the operator k 

lim A = k (28) 

h — >oo 

which corresponds to multiplication on |k| in Fourier space while the operator R for two-dimcntsional flow turns 
into the Hubert transform: 

lim R = H, Hf(x) = -P.V. / ±^-dx'. (29) 

/woo 7T J-oo X — X 

H can be also interpreted as a Fourier transform of i sign(fc). 

If one neglects gravity and surface tension, g = 0, a = 0, then Eqs. ©7 a t leading order over small parameter 
[Vr?|, result in[ll|lllla 

Y t = (30a) 
§ = ^[(^) 2 -(V*) 2 ]- (30b) 

Remarkable feature of Eqs. (|30a|l . (|30b|) is that the second Eq. i|30b[) does not depend on 7/ thus one can first solve 
(I30b|> and then find 77 from Eq. I|30a|) ■ Substitution II = * + iR\fr into Eq. I|30b|) results in the complex Hopf Eq. (gj 
for two-dimensional flow ^3] which is completely integrable. 

Both Eqs. (|30b|l and Q are ill-posed because they have short wavelength instability which is determined as follows: 
we can analyze Eq. (|30b|) and take ^ in the form 

* = * + hie l * r+l/t + c.c), (31) 



where tyo(r,t) is a solution of Eq. (|30b|) . ^1 is the amplitude of small perturbation, and c.c. means complex 
conjugation. Then, in the limit |k| — > 00, if? evolves very slow in space compare to e lk ' r+y * and we get the dispersion 
relation for small perturbations: 

v = A k A-$> - ik- Wo (32) 

which describes instability for Re(v) — A^A^q > 0. For general initial condition such instability region always exists. 
The instability growth rate, Re(v) grows as |k| increases. 

SHORT WAVELENGTH STABILITY ANALYSIS OF THE FOURTH-ORDER HAMILTONIAN 

To study linear stability of the Hamiltonian system in respect to short wavelength perturbations one can set 

»7k = ?7ok + <5f?k, 

*k = *ok + <5*k, (33) 

where 770k, ^ok are solutions of Eqs. Q), (JfJJl and 6r) k , <5-0k are short wavelength perturbations localized around wave 
vector k, |k| 3> q, q is a typical wavenumber for 770k, ^ok- 

If we take into account contribution to the fourth-order Hamiltonan up to second order in amplitude of perturbations 
(5%, i^k we get the following general form of the perturbed Hamiltonian: 

SH = ^J i k |^k| 2 rfk+i J B k \S Vk \ 2 dk + J(F k + l G k )S^ k S v ^dk, 

A k = A_ k , B k = -B_k, Fk = F_k, Gk = — G_k, (34) 
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where A k , _B_k, Fie, G k are real and depend on r/ok, ^ok- Here we disregard linear contribution to SHq because it 
has no effect on linear stability analysis. 

It follows from Eqs. 12U|) . i|34[) that equations of motion take the following form: 

= ik<5*k + (F k - iG k ) Vw , 

at 

= -B k (577 k - (F k + zG k )* k . (35) 
An assumption of exponential dependence on time, 

Sr) k ~ e Vk \ SV k ~ e^*, (36) 

gives a dispersion relation 



v k = -iG k ±\/F 2 - ^kBk (37) 
which describes instability provided F£ — A k B k > 0. 



ILL-POSEDNESS OF THE FOURTH-ORDER HAMILTONIAN 



Consider now a general case of nonzero g and a and take into account all terms in the Hamiltonian up to forth 
order, i.e. consider full Eqs. At the leading order over steepness and wavenumber k we obtain: 

i k = A k + (k 2 - 4)770 - A k {k 2 - A 2 )r, 2 + <3(fc9 2 ), 
F k = B k + A k (A^ ) 2 + 0(k°v 2 l^), 
F k = -A k AV - (k 2 - A 2 k )(A^ )vo + A k (A[ m A^ ] + tj V 2 ^o) 

+O(k Q v a l 1 ), 
G k = k-V*o + O(fcu e), (38) 

where 770 = ^ / r] k dk, ^ Q = ^ J <]/ k cik, and the steepness is defined as 6 ~ |V?7o|. We introduced here the typical 
value of fluid velocity, vq ~ iV^ol and the typical scale, Iq, of variation of vq and rjo : G ~ T)o/Iq , |Vi>o| ~ vq/Iq- 

Eqs. P7|l. give instability growth rate. We consider particular cases. If <r ^ then in the limit |k| — * 00, we 
have 

v k = ±iVa¥, (39) 

i.e. instability is absent. In derivation of Eq. I|39|l we used exponential smallness of expression k 2 — A\ = 
k 2 /(coshkh) 2 ~ Ak 2 exp(— kh) <C 4fc 2 because limit |k| — > 00 implies fcft 3> 1. Thus finite <r makes problem 
©, (EU - © well-posed. 

Note that for finite depth A k < k we could still have instability at finite range of wavenumbers kh <~ 1. In that 
case fc?7o ~ r]o/h <SC 1 because a typical variation of surface elevation, 770, should be small to allow weak nonlinearity 
approximation used throughout this Article. Because krjo <C 1, Eqs. (|37() . (|38() are reduced to 

v k = -ik- Vy Q ±A 1 k /2 [-B k + 2(k 2 -A 2 )(Ay Q ) 2 r ]Q ] 1/2 , kh ~ 1, (40) 
which gives instability provided 

F fc < 2(k 2 - A 2 )(A^ ) 2 Vo, kh ~ 1. (41) 

E.g. instability occurs for g = a = 0: 

i/k = -ik ■ V* ± fe(i* )%^ fe2 " 4*' ( 42 ) 
We can estimate inequality (|41|l as 



(43) 
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where v ~ |V^o| is the typical velocity of fluid. 

It follows from (|43() that instability occurs for large values of vq (because 770 is small). If gravity dominates, 
g > a /h 2 , then ((43(1 gives gh/v 2 , < i] /h but weak nonlinearity approximation implies that 770//1 <C 1 which indicates 
that the kinetic energy strongly exceed the potential energy. Fluid has enough kinetic energy to easily move upward 
at distance ~ h. As a result, at later stage of evolution weak nonlinearity approximation is violated and surface is 
strongly perturbed at scales ~ h. 

If capillarity dominates, g < a/h 2 , inequality l|43() gives cr/(v§?7o) < 1 and the kinetic energy again strongly exceed 
the potential energy. Assume now that, because of instability for kh ~ 1, at later time of evolution the potential 
energy will be of the same order as the kinetic energy, namely, rj v 2 ~ cr9 2 , O ~ \Vt)\. Then a/(v 2 rj ) < 1 results 
in inequality > 1 which again violates weak nonlinearity approximation. Thus for arbitrary relations between g 
and a/h 2 , and for kh ~ 1, the instability is possible for strong enough velocity of fluid and this instability results in 
violation of weak nonlinearity approximation in course of fluid evolution. In that sense there is no surprise that for 
large velocity there is an instability for kh ~ I. This instability is purely physical which leaves problem well-posed. 

Outside capillary scale we can set a = and get from Eqs. 1(37(1 . ((38(1 that zero capillarity makes Eqs. 125(1 ill-posed 
for k — » 00 : 

^ = -ik • Wo ± 2 1 / 2 (i* ) 1 / 2 fc(i[r/ i*o] + ^oV 2 *o) 1/2 

ik- W ±kv ® 1/2 - (44) 

An expression A^ (A^rjoA^o] + ? 7oV 2 5'o^ in Eq. (|44|l is not sign-definite which results in instability of the system 

(|20ll , (1211) — 1(23(1 ■ It means that the fourth order Hamiltonian does not prevent short- wavelength instability but 
makes instability weaker by the small factor O 1 / 2 compare with instability of the third-order Hamiltonian (compare 
Eqs. (|44|) and l|32|)l. Instability l|44|l has been observed numerically 01 . We conclude that full fourth order system 
PU)|. l|2T|) - is ill-posed for zero capillarity, a = 0. 

Ill-posedness of Eqs. I|2(J|) . I(21|) — l|23|l can be also interpreted as violation of perturbation expansion J2J f° r k ^ 00. 
Namely, short wavelength contribution to the quadratic Hamiltonian (|34() , 1(38(1 is not small compare with the other 
terms in the Hamiltonian 1)21(1 — ((23(1 provided krjo > 1. 

Ill-posedness makes Eqs. ((20(1 . 121(1 — ((23(1 (or, equivalently, Eqs. ((25(1 ) difficult for simulations. There a few 
ways to cope with that problem. One way is to resolve all scales down to capillary scales which is extremely costly 
numerically. E.g., if we want to study water waves in gravitation region (scale of meters and larger), we would have 
to simultaneously resolve capillary scale ~ 1cm. Other way is to introduce artificial damping for short wavelengths, 
i.e. to replace Eqs. (|2"U|) by 

dVk 5H , , dl^ SH 

-f = ^+7x(% k , = - (45) 

where functions 7i(fc), ^(k) are zero for small and intermediate values of k but they tend to —00 for k — > 00. Also it is 
possible to introduce finite viscosity of the fluid. However in that case we would have to resolve very small scales and, 
in addition, the Hamiltonian is not conserved for finite viscosity so that we can not use the Hamiltonian formalism. 

In this Paper we use another way which is to completely remove short wavelength instabilities and make problem 
well-posed by appropriate canonical transform from variables 77, \1/ to new canonical variables £, R. 

CANONICAL TRANSFORM 

Canonical transform from variables ^,77 to new variables R, £ is determined by a generating functional S : 

S = J Rki]- k dk+ -|- J A 3 ij kl rj k2 R k3 5(ki +k 2 + k 3 )dkidk 2 dk 3 

+ ., n , 9 [ Vk 1 .k,.k 3 .k 4 -Rk 1 »7k,?7k3'7k4^(ki + k 2 + k 3 + k 4 )dkicik2dk3<ik4, (46) 
4(2?r)^ J 
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SS 



^2 J ^k 1 ,k 2 ,k 3 ,-k-Rki'7k 2 ??k 3 ^(ki + k 2 + k 3 - k)dkidk 2 dk 3 , 
£k = ^ — = Vk + J ^fc^ki^d^ki + k 2 - k)dkidk 2 
J ^-k,k 2 ,k3,k4 7 7k2 f 7k 3 ? 7k4 ,5 ( k 2 + k 3 + k 4 - k)dk 2 dk 3 dk 4 , 



(47a) 



(47b) 



where Vki,k 2 ,k 3 ,k 4 is the symmetric function of k 2 ,k 3 ,k4. This is the most general form of canonical transform up to 
terms of the fourth order. The only condition which we use here is that S is chosen to be linear functional of R to 
preserve the quadratic dependence of the Hamiltonian on canonical momentum R. 

The quantity r] can be found from Eq. I|47b|) as the functional of £ by iterations (here and below we take into 
account only corrections up to the fourth order in the Hamiltonian): 



Vk 



If 1 
£k - ^ J A-£ kl £k 2 £(ki + k 2 - k)dk lC ?k 2 + — 



A k A 



2V- 



k,ki,k 2 ,k 3 



(2tt)2 

Ck^k^ka^ki + k 2 + k 3 - k)dkidk 2 dk 3 , 



(48) 



Eqs. l|47a|) . give: 



*k = i?k + ^ J A 1 i? kl ^ k2 <5(k 1 + k 2 - k)dkidk a 



8(2^) 2 



- ^llA 2+3 + GVk^ka.ka ,-k 



i?k!^k 2 ^k 3 



xtf(ki + k 2 + k 3 - k)dkidk 2 dk 3 



Using Eqs. J23>, (E3>, {HE) , BE) we get: 



Hi = — 



1 

4tt 



H, = 



Hq = IJ { A kW + B k \^\ 2 }dk, 

(kx • k 2 )i? kl i? k2 - -{AxBt + A 2 B 2 + A 3 B 3 )^ k2 
6 J 

x£ k3 <5(ki + k 2 + k 3 )dkidk 2 dk 3 , 

(kx • k 2 )(A 1+2 - A 1 - A 2 ) - k\A 2 - k\A x 



(49) 



(50) 



(51) 



8(2^) 2 



+ -AiA 2 [v4i +3 + A 2+3 + A 1+i + A 2+i ] + 3[AiT/ k2 , k3 , k4jkl 
-^2V kl!k 3 !k4 , k2 ]|i? kl i? k2 ^ k3 ^ k4 (5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 dk 4 

+ 8(2^)2 / { ~ CT(kl ' k2)(k3 ' k4) + i j4 ?+2 Bl + 2 + A 3B 3 A 1+2 

- 2Bi ,k 2 : k 3 ,k 4 } £ki 6c 2 6c 3 £k 4 

x<5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 (ik4, 

Bj = B kj , B J+ i = Sk 3 +k r 



(52) 



Canonical transform conserves the Hamiltonian structure so the dynamical equations in new variables R, £ are 
given by: 



<9£ _ £ff di? 5Jf 



(53) 



8 



FROM COMPLEX TO REAL HOPF EQUATION 



We choose the cubic term of the generating functional S in such a way to remove linear instability at leading order. 
Similar to Eqs. j3()ajl . I|30b(l . we get from Eqs. fKHjl . I|51|l. (|53|l at leading order of small parameter |V£| : 



(54a) 
(54b) 



Thus instead of the complex Hopf Eq. (or Eq. (|30bjl ) we got the real Burgers Eq. I)54bl) for new canonical variable 
R. It is important that the real Burgers Eq. is well-posed. 

Additional advantage of Eq. I|54b|) is that it can be integrated by the the method of characteristics not only in two 
dimensions as Eq. but for three dimensional flow also. 



REMOVAL OF INSTABILITY FROM FOURTH ORDER TERM 



Next step is to remove instability from the fourth order terms in the Hamiltonian Q52J1 by a proper choice of matrix 
element V . We can take Vki,k 2 ,k 3 ,k 4 in the following form: 



Vk!,k 2 ,k3,k 4 = a.\k\ + a 2 Ai(A 2+3 + A 2+i + ^3+4), 
where a.\, a 2 are the real constants. The Eqs. (|50|l . (|51l) . 1521) take the following form: 

1 



Hn = 



1 



■ / {A k \R*\ 2 + B k \&\ 2 } 



- (kx • k 2 )i? kl .R k2 - -{A 1 B 1 + A 2 B 2 + A 3 B 3 )^ k2 
x^ k3 <5(ki + k 2 + k 3 )dkidk 2 dk 3 



Ho 



(k 1 -k 2 )(A 1+2 -A 1 -A 2 ) 
+ (3ai - l){k\A 2 + k\A x ) 



+ [ 4 + 3a 2] AxA 2 [A 1+3 + A 2+3 + A 1+i + A 2+i ] 

+3a 2 A 1 A 2 (A 1+2 + A 3+4 )} 
xi? kl i? k2 ^k 3 ^ k4 (5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 dk4 

"8(hy J { ~ CT(kl ' k ' )(k3 ' k4) + I A i+ 2Bl + 2 + A 3S 3 A 1+2 
-2Bi [ ai kl + a 2 A 1 (A 2 +3 + A 2+4 + A 3+4 )] 

x<5(ki + k 2 + k 3 + k4)dkidk 2 dk 3 (ik4, 



Ba = Bi, 



B 



B) 



kj+k, 



(55) 



(56) 



(57) 



(58) 



9 



The dynamical equations, as follows from (|53() . (|56() — (|58|l . are 



^Vi?- VAR- -(1 - 3ai)£(Afl)AR- (- +3a 2 )(AR)A(£AR) 



a\ A ^ 3 3a 2 



--^(AB^-^Ai^ 



<K -AR-V- [(VE)e]+|V 



-iv-(C 2 iVi?) + ^^iA(C 2 ii?) 



4 
1 

~ 4 

1 — 3ai 



AE(£A£ 2 ), 
V • i(£ 2 Vi?) 
A(£ 2 Ai?) 



1 



3a 2 )A £i(£AR) 



-a 2 A 



(i^ 2 )(ii?) 



(59) 



(60) 



where 5 = «? - <rA, A = ^ + ^. 

To study linear stability of the Hamiltonian system in new variable in respect to short wavelength perturbations 
one can set, similar to Eq. (|33J) . variables £, R in the following form 



£k = £ok + <5£k, 
^?k = -R()k + 5Rk, 



(61) 



with an assumption of an exponential dependence on time 

<*& ~ e» kt , SR k 



(62) 



Here £ok, -Rok are solutions of Eqs. (|59|l . (|60|l and (5£ k , <5^?k are short wavelength perturbations localized around wave 
vector k, |k| 3> q, q is a typical wavenumber for £ok, ^?ok- 

We get, similar to Eqs. 134|l . I|35(l . (I37|l . for the perturbed Hamiltonian, 



A k |«? k | 2 dk 



Bk|^k| 2 dk + / (F k + iG k )«? k ^_ k dk, 



A_ 



£>- k , F k — -F_ k , G k 



-G_ 



(63) 



the following expressions: 



A k = A k + k 2 £ + A k [^ ai k 2 + (i + 3a 2 )A 2 ]£ 2 + O(fc 2 6£o), 
B k = B fc - A k B k Z + (- + 3a 2 )A k (AR ) 2 

+ZB k [(\-a 2 )Al~\aik 2 M 
+O(gke£o) + 0(ak 3 e^) + O(fc% V), 
F k = [ - ifc 2 (l - 3a!) + (i + 6a 2 )A 2 ] (ARo)£o + 0(^6), 

G k = k- Vi? o + O(fcw o 0), 



(64) 



where £,0 = ^1 Ck^k, Ro = ^ I Rkdk, and 9 ~ |V£o| is the steepness. Similar to Eq. (|38|) . we introduced here the 
typical value of fluid velocity, vo ~ |Vi?o| and the typical scale, /o of variation of Do and £o : © ~ £oAo> |Vwo| ~ vq/Iq. 

Eqs. I(37jl . (164(1 give instability growth rate. Our purpose is to make these Eqs. well-posed for zero capillarity so 
that we assume a = and consider limit k — > oo which means that A k — k. It is convenient to rewrite Eqs. (|37|l . (I64|l 
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in dimensionless form as follows: 



A k = 1 + A + [~ai + i + 3a 2 ] A 2 + 0(A6), 
^k = 1 - A + (~ + 'ia 2 )p + 3[i - a 2 - ^] A 2 + O(A0) + 0(A"Ve), 

f£ = [ - \ + \ax + 6a 2 ] 2 A 2 p + 0(Ap6), 



-2 



H - ^k^k, (65) 



where A k = A k /k, _B k = B^/g, F£ — F^/(gk), = (z^k + iG k ) 2 /(gk). The system (IfiBT is described by the two 
independent dimensionless parameters A = fc£o and p = k(ARo) 2 / g which reflects the freedom of choice of an initial 
surface elevation and an initial velocity. Condition of applicability of Eqs. (|64() is Mq <C 1, which gives |A| ^ in 
dimensionless variables p, A. Parameter p can take any nonnegative value because it depends on the fluid velocity 
which can be arbitrary. We want to choose a\ and a 2 to ensure that Eqs. J^SJ, (EEJ, ijBTf . (p>H)l are well-posed and 
stable, which means that v 2 ^ < 0, for any value of A and any nonnegative value of p. 

First step is to analyze the system (I65f) in the limit O — > which means that we first neglect 0(. . . O) terms in l|65|l . 

Assume that g ^ then the necessary condition for v 2 ^ < is to have Ak-Bk > which means that either Au > and 

_/3 k > 0, or A k < and -B k < 0. It easy to show that the second case can not be realized for Eqs. (|65[) so we consider 

the first case of positive A k and _B k . Inequality A k > gives 0% > —1/3 and inequality B k \ p= o > gives (3i < 1/6 
which together result in 

-^ </?1< lT' (66) 

Provided (|66|) is satisfied, the sufficient condition for absence of instability, < 0, is to have term oc p in v? to be 
negative for any A, which means that 

[ - \ + l<*i + 6 «2] 2 A 2 3<* 2 )i k < 0. (67) 

This inequality is satisfied for any A provided 

Pi+P2<^, /3 2 = J(-l + 3a 1 + 18« 2 ). (68) 

Thus v 2 . < for any A and p > provided inequalities (|66|l and (|68|l hold. It corresponds in (a>i, a 2 ) plane to the 
inner part of the ellipse defined by H68J) and bounded by two parallel lines defined by l|66(l (see the filled area in Figure 
1). The center of the ellipse is located at a% = a 2 = ^ (point A in Figure 1). So the choice of a.\ and a 2 is not 
unique for g ^ in the limit — > and is determined by (|66|) , (|68|) . 
For 5 = Eqs. I|65|l are reduced to 

Ik = 1 + A + [^ai + J + 3a 2 ] A 2 + 0(A6), 
~vl = [- J + ^ai+6a 2 ] 2 A 2 - [j + 3a 2 + 0(\~ 1 Q)]A k + 0(A6), 

^k = & ^k + *Gk) 2 - (69) 

A 2 (AR ) 2 

It follows from Eq. Ij69(l that v 2 ^ < in the limit O — > provided condition (|68(1 is satisfied, which corresponds in 
(ai, a 2 ) plane to the inner part of the ellipse defined by (|68|l (see Figure 2) in contrast with the case of nonzero 
gravity. 

For fc-toowe get from (|64|) for g ^ 0: 



^ k = -ik ■ VR ± ik^g 1 ^ 2 + J + 3a 2 ] 1/2 



x[i- a2 -f] 1/2 Co 2 (70) 
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and 

z/k = -ik ■ Vi? ± ik 2 2~ 3/2 3 
x [ - 2a\ + 4a 2 (l - 6a 2 ) + ai(l - 12a 2 )] 1/2 £ ii?o (71) 

for 5 = 0. 

Note that for parameters satisfying inequalities (|66[1 . (|68|1 . the real part of is zero even for kh ~ 1, where ^k takes 
the following form: 

^ k = zk • Vi? ± *^I /2 [g + ( j + 3a 2 ) A fe (ii? ) 2 ] 1/2 , (72) 

so that in new canonical variables l|55[l . ()66|l . ()68|l the instability is absent even for intermediate values of A: ~ l/h. 
Another remark is that these variables leave problem well-posed for a ^ also but that case is not so interesting 
because Eqs. (|53|) . I|56l) . J57J, (|58|l well-posed even in original variables r\ 1 ip for cr ^ 0. 

Now we make the second step and assume that is small but nonzero in Eqs. I|65[) and i|69|) . Terms 0(. . . 0) in 
Eqs. i|65|) and i|69|) are not sign-definite and their values depends on horizontal coordinates x, y and time according 
to dynamical Eqs. I|53|) . i|56[) . (|57|l . (|58l) . Generally these terms result in shrinking of the area of stability, vfc < in 
(ai, 012) plane. Figures 3 and 4 show shrinking of the stability area for the particular choice of terms 0(. . . 0) for 
g 7^ 0. Each curve in Figure 4 corresponds to the stability boundary, max^ = 0, for the particular value of 0. The 

additional requirement is that 

^0(2 + hS) <a 2 + Y<\ + ^©(2 " 626), (73) 

which is a generalization of l|66|> for nonzero 9. Here we assume O(A0) = 61 A8 and O(A0) = fe 2 A0 in rhs of Eqs. 

(16511 for Ak and B^, respectively. Calculating curves in Figures 3 and 4 we set b\ — b 2 — — 1 as a typical example. 
Condition (|73[1 result in additional cutting of curves max^ = in Figures 3 and 4 for = 0, 0.01, 0.02, 0.025 and 

A, p 

0.05. The system (|65|1 is stable inside each curve in Figure 4 for given 0. For < 0.1, the width of region between 
solid curve (0 = 0) and curves with 0^0 scales approximately as 0. For > 0.1 the region of stability quickly 
shrinks to zero as approaches ~ 0.1389. Note that these numerical values are non-universal and depend on the 
numerical coefficient in 0(. . . 0) terms. In a similar way, Figures 5,6 show shrinking of the stability area for zero 
gravity case. 

Our objective is to find parameters a\, a 2 corresponding to stability, v\ < 0, with the largest possible 0. For the 
case g 7^ this is achieved if the maximum of ^|e=o, as a function of A, p, is not only negative but minimum as a 
function of a\ and a 2 , i.e. we want to find min maxi^ „ . This ensures that the system (|65J) is the most stable 

01,02 A, p IW— U 

system for 0^0 or, in other words, the system is the most rigid one. Because p > we have to set p = to find 
maxi^. Then we obtain that min max^ = v&\ . = — 1. This minimum is attained provided 3ai + 6a 2 = —1/2 

A,p ai,a 2 A, p lp— U, A— U 

and i(l — y5) < cti < 1/3. We also want to have the most stable system for k — > 00 (p — » 00 and A — > 00). This is 
achieved provided the coefficients in front of the leading order terms A 4 and X 2 p in v\ are minimums. The coefficient 
for A 4 is already minimum from condition 3ai + 6a 2 = 1/2, while the coefficient for X 2 p is —3/16 + 9a 2 /A, i.e. we 
have minimum for 

ax =0, a 2 = 1/12, g ? 0, (74) 

which corresponds to point B in Figure 1. This choice of parameters a.\ and a 2 is optimal to keep the system l|65f) 
stable for the largest possible 0, i.e. for the largest possible nonlinearity. 

In a similar way, for zero gravity, g — 0, the system l|69|l is the most stable provided we find a±, a 2 which correspond 
to min max^L for Eqs. I|69|l . Maximum rnaxi/^ „ is attained for 

01,02 A IW— U ^ 1 W— U 

A = ^ 1+ ^ 2a2 ^ [2a 2 + 4a 2 (6a 2 - 1) + a 1 (12a 2 - 1)] ~* , 
9 

and min max^L = I (—2 + y3) is attained provided 



kle=o 3 



"i=-2+3- 1/2 , a 2 = i -3 3 / 2 , g = Q, (75) 
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which corresponds to point E in Figure 2. This choice of parameters ot\ and a 2 is optimal to keep the system 16511 
stable for the largest possible 9, i.e. for the largest possible nonlinearity. 

Thus we can choose a±, a 2 from the conditions l(66|) and ((68(1 to make Eqs. ((53(1 . (|56|) — (|58|) (or, equivalently, Eqs. 
(|59H . (|60|l ) well-posed for any value of cr, g and arbitrary depth of fluid. To find dynamics of free surface, one can solve 
Eqs. for R, £ using Eqs. (|56|) — (|58l) and conditions (|66|) . I(68|) . This is the main result of this Article. To recover 
physical variables ^, n from given i?, £ one can use Eqs. l(4*S|) . l(^$|) . ljS3|l . 

Now we can return to the comment in Section 5 about interpretation of ill-posedness of Eqs. (|20|l . (|21|1 — 1(23(1 as 
violation of perturbation expansion (J2J) for krjo > 1. For the new canonical variables £, R perturbation expansion 
is still formally violated for k£o > 1 because contribution from the quadratic Hamiltonian ((63(1 . ((64(1 is not small 
compare with other terms in the Hamiltonian ((56(1 — 1(58(1 . However this violation does cause any problem because 
there is no short wavelength instability in the new canonical variables and the system l(53|l , l(56|l — (|58|) is well-posed. 
In other words, the new canonical variables £, R provide purely physical way to regularize short wavelengths without 
introduction any artificial viscosity. 

As follows from Eqs. 166(1 . 1(68|l . the new canonical variables £, R are not uniquely determined from the condition 
of well-posedness of the dynamical Eqs. (JHHJl , (J23) — because parameters (ai, a 2 ) can take any valued from 
filled area in Figures 1 and 2. However the choice of (ai, a 2 ) is unique provided we additionally require the system 
((5311 . (|56l) — H58fl to be free of short wavelength instability for the largest possible slopes 9, i.e. for the largest possible 
nonlinearity. This gives the conditions (|74|1 for g 7^ and 1(75(1 for g = 0. We refer to the variables £, R, defined in 
Eqs. 1(48(1 . 1(49(1 . 155 . [I. ((74(1 and 1(75(1 . as the optimal canonical variables. For some extent similar results were obtained 
by Dyachenko 18] for particular case of two-dimensional flow. We conjecture that the optimal canonical variables, 
which allow well-posedness of the dynamical Eqs., exist in all orders of nonlinearity However additional research 
necessary to decide if the optimal canonical variables exist and unique in higher (fifth etc.) order of nonlinearity. We 
also conjecture that the optimal canonical variables £, R would allow simulation with higher steepness compare with 
standard variables n. 



There are a number of important special cases of the optimal canonical variables. Here we use ((74(1 and 175(1 . We 
give here expression for the Hamiltonian only. The dynamical Eqs. can be obtained either from ((53(1 or directly from 



SPECIAL CASES 



Eqs. (ED, 



Deep water limit 



For g ^ 0, h — > 00, Au = k and Eqs. (|55|l 



1(58(1 take the form 




(76) 



(77) 



+ 



8(2^) 2 



1 



/{ 



xi? kl i? k2 £ k3 £ k4 <5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 (ik4 
<r(ki • k 2 )(k 3 • k 4 ) + -|ki + k 2 | 2 B 1+2 + /c 3 £ 3 |ki + k 2| 




xCkiCks&y&^ki + k 2 + k 3 + k 4 )dk 1 dk 2 dk 3 dk 4 . 



(78) 
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Zero gravity and capillarity g = a = 



H ° = \J k\R*\ 2 dk, (79) 



Hi = / (ki • k 2 )i? kl i? k2 e k3 ^(ki + k 2 + k 3 )dk 1( ik 2( ik 3 , (80) 



47T 



^ 2= 8(ty J {( k i- k 2)(|ki+k 2 |-fc! 



k 2 ) 



+ (-^ + V3)(fc 2 fc 2 + fc 2 2 fci) 

+ (1 - 3 5/2 )fc 1 fc 2 [|k 1 + kg | + |k 2 + kg I + |ki + k 4 | + |k 2 + k 4 |] 

+ (^-3 5 / 2 )fc 1 fc 2 [|k 1 +k 2 | + |k 3 +k 4 |]} 
xi? kl i? k2 ^ k3 ^ k4< 5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 <ik 4 . (81) 



Zero gravity, g — 0, and nonzero capillarity a / 

Ho = \J {k\R*\ 2 + o-k 2 \&\ 2 }dk, (82) 

Hi = J [ - (ki • k 2 )i? kl i? k2 - + k\ + fc|)a^k 2 ] 

x^ k3 5(ki + k 2 + k 3 )dkidk 2 dk 3 , (83) 

H2 = s(ky / ( (kl ' k2)(|kl + ka| ~ kl k2) 

+{-\ + VTi){k 2 k 2 + fc 2 2 fci) 

+ (1 - 3 5/2 )fc!fc 2 [|ki + k 3 | + |k 2 + kg I + |ki + k 4 | + |k 2 + k 4 |] 

+ (^-3 5 / 2 )fc 1 fc 2 [|k 1 +k 2 | + |k 3 +k 4 |]} 
xi? kl i? k2 ^ k3 ^ k4 (5(ki + k 2 + k 3 + k 4 )dkirfk 2 dk 3 rfk 4 

+ &kf J { ~ (kl ■ k2)(k3 ■ k4) + I |kl + k2 ' 4 + * + kz| 

-kf[(-l + Z- 1 ' 2 )^ + (1 - 3 3 / 2 2)(|k 2 + kg | + |k 2 + k 4 | + |k 3 + k 4 |)] } 

x£ kl £k 2 £ k3 £ k4 <5(ki + k 2 + k 3 + k 4 )rfkidk 2 rfk 3 <ik 4 . 



(84) 
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Nonzero gravity, g 7^ 0, and zero capillarity a = 



Hi 



4tt 



Ha = \j {^k| 2 +5l&| 2 }dk, 

(ki • k 2 )i? kl -Rk 2 - ^(h + fc 2 + k 3 )£ kl £k 2 
6 

x£k 3 <5(ki + k 2 + k 3 )dkidk 2 <ik 3 , 



ff2 ^8(^F/{ (kl - k2)(|kl + k2 



I - ki - fc 2 ) - (fc 2 fc 2 + fe 2 fc i) 



+-fcife 2 [|ki + kg I + |k 2 + kg I + |ki + k 4 | + |k 2 + k 4 |] 

+ifc 1 fc 2 [|k 1 +k 2 | + |k 3 + k 4 |]} 
xi? kl i? k2 ^ k3 ^ k4 (5(k 1 + k 2 + k 3 + k 4 )dkidk 2 dk 3 cfl<4 



8(2^) : 

-^ x (|k 2 + ka| + |k 2 + k 4 | + |k 3 + k 4 |)} 
x£ki£k 2 £k 3 £k 4 <5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 dk 4 . 



(85) 



(86) 



(87) 



Shallow water limit 



Shallow water limit corresponds to kh — > 0. In that limit — > fe 2 /i. Eqs. i|56|l — (|58|) take the following form for 
.9^0: 



Hi = — 

4tt 



H ° = IJ { k2h \Rk\ 2 + Btl^dk, B k =g + ak\ 

- (k x • k 2 )fl kl fl ka - + fc|B 2 + /c 3 2 B 3 )£ kl 6c 2 

x^ k3 <5(ki + k 2 + k 3 )dkidk 2 dk 3 , 
J {2(ki • k 2 ) 2 - + fc 2 fc 2 ) 



8(2tt) : 



+/i 2 -fc 2 fc 2 [|ki + k 3 | 2 + |k 2 + k 3 | 2 + |kj + k 4 | 2 + |k 2 + k 4 | 2 ] 

+/i 2 ifc 2 fc 2 2 [|k 1 +k 2 | 2 + |k 3 + k 4 | 2 ]} 
xi? kl i? k2 Ck 3 ^k 4 ^(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 dk 4 

+ 8( 2 V/{-' j(kl - k2)(k3 - k4) 

h 2 



^ |ki + k a |*B 1+ 2 + ft^B 3 |ki+k 2 ^ 

(|k 2 + k 3 | 2 + |k 2 + k 4 | 2 + |k 3 + k 4 | 2 )} 
x£ kl £ k2 £ k3 £ k4 <5(ki + k 2 + k 3 + k 4 )dk!dk 2 cfk 3 (ik 4 . 



6 



(88) 



(89) 



(90) 



Zero gravity and capillarity g = a = 



H = \ J k 2 h\R^\ 2 dk, 
Hi = J (ki • k 2 ) J R kl flk 2 & 3 <K k i + k 2 + k 3 )dkidk 2 dk 3 , 

ff2 = / ^ 2(kl ' k2)2 + + ^ ){klkl + klkl) 

+h 2 (l - 3 5 / 2 )k 2 k% [|ki + k 3 | 2 + |k 2 + k 3 | 2 + |ki + k 4 | 2 + |k 2 + k 4 | 2 ] 

+ /,2 ( 3 _ 3 5/2 )fc 2 fc 2 [|ki + k2| 2 + |kg + k4|2] | 

xi?k 1 i?k 2 Ck 3 Ck 4 <5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 (ik4. 

Zero gravity, g — 0, and nonzero capillarity a / 



° 2 



ffi = ^ / [ - (ki • k 2 )i? klJ R k2 - ^(tf + fe| + fc 3 4 )? ki a 2 ] 



x^ k3 (5(ki + k 2 + k 3 )rfkidk 2 (ik 3 , 



ff 2 



8(2tt) 2 

+fe 2 (l-3 5 / 2 )fc?A ; |[|k 1 +k 3 
+^ 2 ( 



2(ki • k 2 ) 



5 -•••'•2 .2.2 , 



(__ + V3)(fc 2 fc; 



2,3 



|k 2 +k 3 | 2 + |k 1 +k 4 | 2 
-3 5 / 2 )fc 2 fc 2 [|k 1+ k 2 | 2 + 



k^k\ 

-|k 2 + 
k 3 + k 4 



k 4 r 

2 



1} 

X-Rk^k^kaCk^ki + k 2 + k 3 + k 4 )dki<ik 2 dk 3 (ik 4 

+ 8(2^ / { ~ ^ ' k2Kk3 ' k4) + T |kl + k2 ' 6 + ^s^ 1 + k2 ' 2 
-kf [(-1 + 3-V 2 ) + (I - 3 3 / 2 2)/i 2 (|k 2 + k 3 | 2 + |] 



k 4 | 2 



k 4 | 2 ) 



<^k 1 Ck 2 Ck 3 ?k 4 ^(ki + k 2 + k 3 + k 4 )dkidk 2 rfk 3 dk 4 . 
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Nonzero gravity, g 7^ 0, and zero capillarity a = 




(97) 



Hi = ~r / [ - (k x • k 2 )i? kl i? k2 - ^r(fc 2 + fc 2 2 + fc 3 2 K kl & 2 




x^ k3 (5(ki + k 2 + k 3 )dkidk 2 dk 3 , 



(98) 



H2 = ^ / {2(k x • k 2 ) 2 - (fc 2 fc 2 + fcffc 2 ) 

+/ l 2 ifc 2 fc 2 [|ki + k 3 | 2 + |k 2 + k 3 | 2 + |ki + k 4 | 2 + |k 2 + k 4 | 2 ] 
+/ l 2 ifc 2 fc 2 [|k 1 + k 2 | 2 + |k 3 + k 4 | 2 ]} 



xi? kl i? k2 Ck 3 Ck 4 ^(ki + k 2 + k 3 + k 4 )dki(ik 2 dk 3 dk4 




x£ki£k 2 £k 3 £k 4 <5(ki + k 2 + k 3 + k 4 )dkidk 2 dk 3 (ik4. 



(99) 



CONCLUSION 



In conclusion, we found the optimal canonical variables for which the water wave problem is well-posed in the 
approximation which keeps terms up to fourth order in the Hamiltonian. The important question remain open if it is 
possible to make water wave equations well-posed by proper choice of canonical transform for higher-order corrections 
(fifth and higher order). We conjecture that such optimal canonical variables exist in all orders of nonlincarity. 
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0.15- 




FIG. 1: Eqs. (I53H . 1561 — 1581 with j / are stable for any (ai, an) inside the filled area. Point A corresponds to the 
centre of ellipse, (ai, a.2) = (|, gg). Point B corresponds to the most stable system 16511 . (at, a.2) = (0, ^). Points C, 
(«i, 0L2) = (5, — Y2), and D, (ai, Q2) = (g, correspond to intersections of the ellipse defined in 168H with the line 
3ai + 602 = 1. lushnikovzakharovfigl.eps 
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0.15- 




FIG. 2: Eqs. 1531 . 1561 — 158H with g — are stable for any (on, a^) inside the filled area. The filled area is bounded by the ellipse 
defined in 1681 . Point A corresponds to the centre of ellipse, (ai, oti) = {\, Point E, (ai, a.2) = (—5 + 3 -1 ' 2 , | — 3 3,/2 ), 
corresponds to the most stable system 1691 . lushnikovzakharovfig2.eps 




FIG. 3: Shrinking of the stable region for small but nonzero O for g 7^ 0. Solid curve corresponds to O = 0. Area inside dotted 
curve corresponds to the stable region for Q = 0.025. Dotted curve is obtained numerically for Eq. (1651 where we set as example 
O(A0) = — X&, 0(X~ 1 pQ) — X~ 1 p<d, O(XpO) — —XpO. Points A, B are defined in Figure 1. lushnikovzakharovfig3.eps 
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6=0.00 
6=0.01 




FIG. 4: A series of curves corresponding to shrinking of the stable region as a function of for j / 0. For O < 0.1, the 
width of region between solid curve (O = 0) and curves with 6/0 scales approximately as G. For > 0.1, the stable region 
quickly shrinks to zero as approaches ~ 0.1389. Similar to Figure 3, all curves are calculated numerically with assumptions 
O(A0) = -A0, O(A _1 p0) = A~Ve>, O(Ap0) = -Ap0. lushnikovzakharovfig4.eps 
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a 



FIG. 5: Shrinking of the stable region for small but nonzero Q for g = 0. Solid curve corresponds to O = 0. Area inside 
dotted curve corresponds to the stable region for Q = 0.025. Dotted curve is obtained numerically for Eq. (1651 where we set as 
example 0(A6) = -A9, 0(A _1 9) = A _1 9. Points A, E are defined in Figure 2. lushnikovzakharovfig5.eps 
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6=0.00 
6=0.01 




FIG. 6: A series of curves corresponding to shrinking of the stable region as a function of for g = 0. For O < 0.1, the 
width of region between solid curve (O = 0) and curves with 6/0 scales approximately as G. For > 0.1, the stable region 
quickly shrinks to zero as approaches ~ 0.1389. Similar to Figure 5, all curves are calculated numerically with assumptions 
O(A0) = -A0, O(A _1 0) = A _1 0. lushnikovzakharovfig6.eps 



